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ABSTRACT 

We  discuss  a  subroutine  package  for  solving  the  eigenproblem  of 
bisymmetric  Toeplitz  matrices.  When  the  eigenvalues  and  eigenvectors  of  a 
bisymmetric  Toeplitz  matrix  are  computed  storage  requirements  and  execution 
time  are  reduced  by  taking  advantage  of  the  special  structure  of  the  matrix. 
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INTRODUCTION 


An  n-tn  order  real  symmetric  matrix  M=(MjR)  for  which 


Mjk  “  Mj+1,  k+1  W 

is  called  a  Oisymrnetric  Toeplitz  matrix.  Bisyrnmetric  Toeplitz  matrices  arise, 
for  example,  as  correlation  matrices  in  tne  study  of  discrete  stationary 
random  processes.  Each  element  of  the  correlation  matrix  is  a  function  of 
j-k;  that  is, 

Mjk  =  Mj-k  =  Mk-j 

(e.g.  =  cos(j-k)c,  where  c  is  a  real  number.) 

In  [1],  using  the  Disymmetry  of  matrix  M,  Goldstein  shows  that  the 
eigenvalues  and  eigenvectors  of  matrix  M  can  be  computed  from  two  real 
symmetric  matrices  one-half  its  size,  if  it  is  of  even  order.  In  this 
memorandum,  we  provide  a  FORTRAN  subroutine  package  that  is  an  implementaiton 
of  this  reduction.  By  using  this  package,  one  reduces  both  computer  time  and 
storage  ordinarily  needed  to  calculate  the  eigenvalues  and  eigenvectors  of 
Disymmetric  Toeplitz  matrices. 
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ALGORITHM 

The  sub-program  package  in  the  Appendix  computes  the  eigenvalues  and 
(optionally)  the  eigenvectors  of  an  even  order,  n,  real  symmetric  matrix  M, 
whose  coefficients  satisfy  Eq.  (1),  from  two  real  symmetric  matrices: 


M  jk(+)  ~  Mjk  +  Mj,n+l-k  =  M~kj(+)  ^ 

(k=l,  •••,  n/2;  j=l,  •  •  • ,  k) 

M  jk(-)  =  Mjk  -  Mjjn+i_k  =  (3) 

one-half  the  size  of  matrix  M,  by  [1,  2]. 

The  matrices  M~(+)  and  M~(-)  are  computed  from  the  first  row  elements  of 
matrix  M.  For  example,  by  Property  (1),  Formula  (2)  reduces  to 

M~jk(+)  =  Ml,k-(j-l)  +  Ml,n-(k+j)+2  (4) 

=  ^  kj(+)  ^=1,  •••»  •••»  ^ 


where  addition  of  the  coefficients  of  the  first  row  of  matrix  M  is  replaced  by 
subtraction  when  computing  M~(-). 
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Forming  M~(+)  first,  the  subroutine  package  calculates  the  eigenvalues 
and  (optionally)  the  eigenvectors  of  M~(+).  Then  M~(-)  is  formed  and  its 
eigenvalues  and  (optionally)  eigenvectors  are  computed. 

The  eigenvalues  and  eigenvectors  of  matrices  M~(+)  and  M~(-)  determine 
the  eigenvalues  and  eigenvectors  of  matrix  M.  If  z  is  an  eigenvector  of  M~(+) 
with  corresponding  eigenvalue  d,  and  J  is  the  identity  matrix  of  order  n/2 
with  its  columns  written  in  reverse  order,  then 

[  z  ]  [  z  ] 

M  [  ]  -  d  [  ].  (5) 

[Jz  ]  [Jz  ]. 

Similarly,  if  z  is  an  eigenvector  of  M~(-)  with  corresponding  eigenvalue  d, 
then 


[  z  ]  [  z  ] 

M  C  ]  =  d  [  ].  (6) 

[-Jz]  [-Jz]. 
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STORAGE  CONSIDERATIONS 


Note  that  in  (5)  and  (6),  since  premultiplying  a  vector  by  matrix  J  just 
reverses  the  order  of  the  vector's  components,  all  the  information  about  the 
eigenvectors  of  matrix  M  is  given  by  the  eigenvectors  of  matrices  M~(+)  and 
M~(-);  therefore,  only  n/2  of  the  components  of  each  eigenvector  of  matrix  M 
have  to  be  stored,  and  a  fortiori,  only  n**2/2  entries  (instead  of  n**2)  are 
required  to  store  all  the  eigenvectors  of  matrix  M. 

Furthermore,  since  the  eigenvalues  and  eigenvectors  of  M~(+)  are  computed 
independently  of  those  for  M~(-),  matrices  M~(+)  and  M~(-)  may  share  the  same 
storage  area  for  a  matrix  of  order  n/2.  Therefore,  the  reduction  requires  a 
maximum  of  (3/4)n**2  +  n  storage  locations  to  accomodate  the  matrices  M~(+), 
M~(-)  and  their  eigenvalues  and  eigenvectors,  instead  of  the  2n**2  +  n 
required  when  the  eigensystem  is  solved  directly  from  matrix  M.  In 
particular,  only  (l/4)n**2  +  n  storage  locations  are  used  to  accomodate  M~(+) 
and  M~(-)  and  their  eigenvalues,  if  only  the  eigenvalues  of  matrix  M  are 
required. 


TIMING 


Computational  results  indicate  that  the  eigenproblem  for  a  large  matrix  M 
can  be  solved  by  computer  at  least  four  times  faster  this  way  than  directly 
from  matrix  M.  For  example,  the  following  table  gives  representative 
execution  times  (in  seconds)  to  compute  the  eigenvalues  of  bisymmetric 
Toeplitz  matrices  of  orders  16,  32,  64,  128,  256  by  the  reduction  SUBROUTINE 
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HATRED  (in  the  Appendix)  and  directly  by  the  IMSL  library  SUBROUTINE  EIGRS  [3] 
on  the  VAX  11/780  in  double  precision  arithmetic. 


Order  (n) 

MATRED  (using  EIGRS) 

EIGRS 

16 

.03 

.07 

32 

.13 

.38 

64 

.73 

2.64 

128 

4.97 

21.00 

256 

38.53 

248.63 

The  following  table  gives  representative  execution  times  (in  seconds)  for  a 
single  precision  version  of  SUBROUTINE  MATRED  and  a  single  precision  version 
of  EIGRS  on  the  UNIVAC  1100/62:* 


Order  (n)  Single  Precision  MATRED  Single  Precision  EIGRS 


16 

.02 

.05 

32 

.10 

.33 

64 

.60 

2.24 

128 

4.08 

16.76 

256 

30.16 

129.04 

*  Since  only  a  single  precision  version  of  EIGRS  is  available  in  the  IMSL 
library  on  the  UNIVAC,  the  UNIVAC  version  of  MATRED  is  single  precision. 
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In  multi-program  computer  environments,  variations  in  the  execution 
time  of  a  program  can  occur  due  to  differences  in  machine  load.  Therefore,  in 
order  to  smooth  out  these  variations,  subroutines  MATRED  and  EIGRS  were  run  a 
number  of  times  on  different  matrices  of  order  n  and  the  average  execution 
time  for  each  order  was  recorded  in  the  tables. 


USER  PROCEDURE 


The  calling  sequence  of  SUBROUTINE  MATRED  (double  precision  version) 
is  defined  in  its  commentary  in  the  Appendix.  The  relocatable  code, 

MATRD.OBJ,  is  in  the  VAX  11/780  directory  [MJG]  on  node  2.  The  UNIVAC  single 
precision  version  of  SUBROUTINE  MATRED  is  in  the  UNIVAC  1100/62  (node  2)  file 
MGOLDSTN*MATRD  and  has  the  element  name  MATRD.  It's  calling  sequence  is  the 
same  as  the  double  precision  version  on  the  VAX,  but  none  of  the  arguments  in 
the  sequence  are  double  precision. 


SUMMARY 


A  subroutine  package  for  computing  the  eigenvalues  and  (optionally) 
the  eigenvectors  of  a  bisymmetric  Toeplitz  matrix  is  available  which  reduces 
storage  requirements  and  execution  time  by  taking  advantage  of  the  special 
structure  of  the  matrix.  In  particular,  if  only  the  eigenvalues  of  the  matrix 
are  required  storage  requirements  for  arrays  and  execution  time  are  reduced  by 
75^  when  the  matrix  is  sufficiently  large. 
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APPENDIX 
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SUBROUTINE  MATRED(IVEC,MR1 ,N,IFULST,MTILFL,IFUL,MTILSY,D,Z,IZ,WK) 

******************************************************************** 

A  subroutine  package  for  the  efficient  solution 

OF  THE  E I GENPROBLEM  OF  REAL  SYMMETRIC  TOEPLITZ  MATRICES 

BY 


MARVIN  GOLDSTEIN  AND  WILLIAM  BABSON 

♦  +  ♦  THIS  SUB-PROGRAM  PACKAGE  COMPUTES  THE  EIGENVALUES  AND  *** 

♦  ♦  ♦ (OPTIONALLY)  THE  EIGENVECTORS  OF  AN  EVEN  ORDER  (n)  REAL  SYM-  ♦*♦ 

♦♦♦METRIC  TOEPLITZ  MATRIX;  ♦  +  ♦ 

M  =  M  (  1  ) 

j  k  j  +  1  ,  k+  1 


=  M 

k  j 

♦  ♦  +  FROM  TWO  REAL  SYMMETRIC  MATRICES; 


M-1+  J 
j  k 


M  +  M  =M~(  +  ) 

jk  j  ,  n  +  1  -  k  kj 


(2) 


(  k=  1  .  ....  n/2 ;  j  =  1  ,  ....  k  ) 


M~  (  -  ) 

j  K 


M  -  M  =  M-  (  -  ) 

jk  j.n+l-k  kj 


(3) 


♦♦♦ONE” HALF  THE  SIZE  OF  MATRIX  M,  BY  |1.  2). 


♦  +  ♦  the  MATRICES  M~ O  )  AND  M~ ( - )  CAN  BE  COMPUTED  FROM  THE  ♦** 
♦♦♦FIRST  ROW  ELEMENTS  OF  MATRIX  M.  FOR  EXAMPLE,  BY  PROPERTY  (1),+** 
♦♦♦FORMULA  (2)  REDUCES  TO  ♦ *  * 


M  ~  ( +  )  =  M  +  M 

jk  1  ,  k  -  j  “  1  )  l,n-(k+j)+2 


(4) 


=  M~  (  +  ) 
k  j 


k  =  1 . n/2;  j  =  1  ,  ....  k 


♦♦♦WHERE  ADDITION  OF  THE  COEFFICIENTS  OF  THE  FIRST  ROW  OF  MATRIX  *  *  * 
♦ ♦ +  M  Iv  REPLACED  BY  SUBTRACTION  WHEN  COMPUTING  M~ f - )  .  *  *  * 


♦♦♦  IF  z  IS  AN  EIGENVECTOR  OF  M-(+)  WITH  CORRESPONDING  EIGEN-*** 

♦♦♦VALUE  d,  AND  J  IS  THE  IDENTITY  MATRIX  OF  ORDER  n/2  WITH  ITS  *♦♦ 

♦♦♦columns  written  IN  reverse  order,  then 

I  z  |  l  Z  I 

M  |  !  =  d  [  1  .  ( 5 ) 

|JZ  I  Uz  I 

♦♦♦SIMILARLY,  IF  z  IS  an  EIGENVECTOR  of  mm-)  WITH  CORRESPONDING  *** 
♦♦♦EIGENVALUE  cJ  ,  then 


M 


1  -  JZ 


l  2  1 

I  I . 

i  -  Jz  ] 


(6) 


*  *  ■*  computational  results  Show 

***A  MATRIX  THAT  SATISFIES  (1)  CAN 


THAT  THE  EIGENPROBLEM  FOR  *** 

BE  SOLVED  BY  COMPUTER  AT  LEAST*-*-* 


♦♦♦FOUR  TIMES  FASTER  THIS  WAY  THAN  BY  SOLVING  THE  EIGENPROBLEM 
♦♦♦DIRECTLY  FROM  M  FOR  LARGE  VALUES  OF  n. 


*  *  * 
*  *  * 
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6500  C  ***  FURTHERMORE,  SINCE  THE  EIGENVALUES  AND  EIGENVECTORS  OF  *** 
6600  C  ***M~-(+)  ARE  COMPUTED  INDEPENDENTLY  OF  THOSE  FOR  M-(-),  ONLY  *** 
6700  C  ***A  MAXIMUM  OF  ( 3 / 4 ) n *  *  2  *  n  STORAGE  LOCATIONS  ARE  USED  TO  AC-  *  *  * 
6800  C  ***COMODATE  THESE  MATRICES,  THEIR  EIGENVALUES  AND  EIGENVECTORS,  *** 


b900  C  ***INSTEAD  OF  2n**2  ♦  n  WHEN  THE  EIGENVALUES  AND  EIGENVECTORS  ARE*** 
7000  C  ***EVALUATED  DIRECTLY  FROM  THE  ORIGINAL  MATRIX  M.  IN  PARTICULAR,  *** 
7100  C  *  *  *  0  N  L  V  (  1 / 4 ) n  *  *  2  +  n  LOCATIONS  ARE  USED  INSTEAD  OF  n**2  +  n  TO  *** 
72  I  ***ACCOMODATE  THESE  MATRICES  AND  THEIR  EIGENVALUES,  IF  ONLY  THE  *** 
7300  C  ** *EIGENVALOES  OF  MATRIX  M  ARE  REQUIRED.  *** 
7400  C 

7500  C  ***  THE  'UB-PROGRAM  PACKAGE  ONSISTS  OF  SUBROUTINES  MATRED,  *** 
7bUu  C  ***MTLSYM,  MTlFUL.  A  SUBROUTINE  THAT  SOLVES  THE  E I GENPROBLEM  *+* 
7700  C  *  *  *  O  F  A  REAL  SYMMETRIC  MATRIX  IS  REQUIRED.  IN  IMPLEMENTING,  SUB-*** 

760  j  c  •‘♦■♦routine  matred,  we  use  the  imsl  library  routine  eigrs  1 3 1  to  *** 

7900  C  *  *  *  SOLVE  THE  E I  GEN  PROBLEMS  OF  M- ( ♦ )  AND  M-(-);  HOWEVER,  ANY  OTHER*** 


8 0 0 u  c  ***routine  That  lomputes  the  eigenvalues  and  ioptionally)  the  ei-*** 

8100  C  ***GENVECTORS  OF  A  REAL  SYMMETRIC  MATRIX  MAY  BE  USED.  *** 

8  2  0  0  C 

8300  C  ***  SINCE  ROUTINES  THAT  SOLVE  THE  EIGENPROBLEM  OF  A  REAL  SYM-  *** 

0  4  0 1  C  **+ METRIC  MATRIX  REQUIRE  THAT  THE  MATRIX  BE  STORED  IN  EITHER  *** 

8500  C  +**SYMMETRI  STORAGE  MODE  (LOWER  TRIANGLE  STORED  AS  A  ONE  DIMEN-  **+ 

8600  C  ***8I0NAL  ARRAY  IN  ROw  ORDER)  OR  FULL  STORAGE  MODE,  TwO  SUBROJ-  *‘* 

8700  C  ***TINES  ARE  PROVIDED  FOR  "rHlS  PURPOSE.  SUBROUTINE  MTLSYM  WILL  *** 

8800  C  **+COMPUTE  M-(*  AND  M~ ( - )  IN  SYMMETRIC  STORAGE  MODE,  WHILE  SUB-  *++ 

8900  C  ***TINE  MTlFUL  WILL  COMPUTE  THEM  IN  FULL  STORAGE  MODE.  *** 


9  000  C 
9  10  0  C 

9200  C  REFERENCES 

9  3  0  0  C 

9400  C  **+l.  GOLDSTEIN  MARVIN,  REDUCTION  OF  THE  EIGENPROBLEM  FOP  HER-  *+* 

9  50  U  C  ***MITIAN  PERSYMMETRIC  MATRICES,  MATH.  COMP.,  V.3B,  JAN.  1974  *  +  * 

9600  C 

9700  C  *  *  *  2 .  GOLDSTEIN  MARVIN,  "FURTHER  DECOMPOSITION  OF  THE  PSEUDOIN-  *** 

9800  C  *+*VERSE  AND  E I GENSYSTEM  OF  A  HERMITIAN  PERSYMMETRIC  MATRIX,"  ACM*** 

9900  C  *  *  *  S  I GNUM  NEWSLETTER,  V.  13,  NO.  1,  MARCH  1978.  *** 

10000  C 

10100  C  *  *  *  3 .  "EIGRS,"  IMSL  LIBRARV  FOR  VAX-  11/  780,  EIGTH  EDITION,  1980.*** 

10200  C  **»«»**»**«****»»»»********«*****«***«*****«*****»*«»***»»+**«*<■***+ 

10300  C 

10400  C  <><>*><*•  <s  <><>'  x  ><><>*.><*<  ^  <><><*  <.*<'  <x>  <><.**.?<>*.><•>  o  <><  ^  <  *  <>  •  ^  ^ 

1  U  500  C 

10600  C  -MNPUT  VARIABLES  FROM  CALLING  PROGRAM: 

10700  C 

10800  C  MR  1  DOUBLE  PRECISION  REAL  VECTOR  OF  AT  LEAST  LENGTH  N  IN  THE 

10900  C  CALLING  PROGRAM  CONTAINING  FIRST  ROW  OF  MATRIX  M. 

11000  C  N  INTEGER  VARIABLE  -  ORDER  OF  MATRIX  M 

11100  C  IFULST  INTEGER  VARIABLE  SWITCH  -  SET  *0  ZERO  IF  M- (  +  )  OR  M- ( - )  IS 

11200  C  STORED  IN  SYMMETRIC  STORAGE  MODE  BY  EIgRS;  SET  TO  ONE  IF 

11300  C  FULL  STORAGE  MODE. 

11400  C  MTILFL  DOUBLE  PRECISION  REAL.  TwO  DIMENSIONAL  SCRATCH  ARRAY  OF  ORDER 
11500  C  IFUL,  WHICH  HOLDS  M - ( * )  OR  M-(-)  IN  FULL  STORAGE  MODE. 

11600  C  IFUL  INTEGER  VARIABLE  -  SET  TO  ONE-HALF  THE  DIMENSION  OF  MR  1 

11700  C  IN  THE  CALLING  PROGRAM  IF  IFULST=1:  SET  TO  I  IF  IFULST=G. 

11800  C  MTILSY  DOUBLE  PRECISION  REAL  SCRATCH  VECTOR  DIMENSIONED  AT  LEAST 

11900  C  ((N*  +  2  ♦  2 *N )  ) / 8  IN  CALLING  PROGRAM  TO  HOLD  M- ( + )  OR 

12000  C  M-(-)  IN  SYMMETRIC  STORAGE  MODE  IF  IFULST=0;  IF  IFULST=1 

12100  C  THIS  VECTOR  SHOULD  BE  DIMENSIONED  LENGTH  ONE  IN  THE  CAL" 

12200  C  LING  PROGRAM. 

12300  C  I VEC  INTEGER  VARIABLE.  IF  IVEC  IS  ZERO,  THEN  ONLY  EIGENVALUES 

12400  C  WILL  BE  COMPUTED;  IF  ONE,  THEN  EIGENVECTORS  ALSO. 

12500  C  IZ  INTEGER  VARIABLE  -  ROW  DIMENSION  OF  Z  ARRAY  IN  CALLING  PRO- 

12600  C  PROGRAM  THAT  HOLDS  EIGENVECTORS  OF  M- ( + )  AND  M~ ( - )  .  IZ  IS 

12700  C  ONE-HALF  THE  DIMENSION  OF  MR  1  IN  THE  CALLING  PROGRAM  IF  IVEC  = 

12800  C  1;  OTHERWISE,  IZ  IS  ONE. 

12900  C  WK  DOUBLE  PRECISION  WORKSPACE  VECTOR  OF  AT  LEAST  LENGTH  N. 

13000  C 

13100  C  -  >  OUT  PUT  TO  CALLING  PROGRAM: 

13200  C 

13300  C  D  DOUBLE  PRECISION  REAL  VECTOR  OF  AT  LEAST  LENGTH  N 

13400  C  IN  THE  CALLING  PROGRAM,  WHICH  HOLDS  THE  EIGENVALUES  OF 

1 3500  C  MATRIX  M. 

13600  C  Z  DOUBLE  PRECISION  MATRIX  OF  IZ  ROWS  AND  2*IZ  COLUMNS  IN 

13700  C  THE  CALLING  PROGRAM.  ONTROL  IS  RETURNED  TO  THE  INVOKING 

13800  C  PROGRAM  WITH  THE  EIGENVECTORS  OF  M~(  +  )  IN  COLUMNS  1 . N/2 

1 390  u  C  OF  MATRIX  Z,  AND  THE  EIGENVECTORS  OF  M-(-)  IN  COLUMNS 

14000  C  N/2+1 . N.  WITH  CORRESPONDING  EIGENVALUES  IN  VECTOR  D. 

14100  C 

1  4  2  0  u  C  <><><><><><><>*><>«.**»*•*«•.'«>*•  ><•><>»  >*  ><  »<  •«..•<.><><><><><><.><><>«  *  <  »  ‘ 
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14  3  0  0 
14400 
145  J  ) 
1  460U 
14700 
1  4  3  U  U 
1  4900 
1  5000 

15  100 
1  5  2  0  0 

15  30  0 
1  5400 
1  5500 
1  5b00 
15700 
’  5600 
1  5900 
i  6000 

16  100 
1  6  200 
1  6300 
1  6400 
1  6500 
1  to  b  0  0 
1  6  7  00 
1  6800 
1  6900 
1  7000 

17  100 
1  7200 
17300 
1  7400 
1  7500 
1  7  600 
1  7  700 
1  7800 
1  7900 
18000 

18  100 
1  8200 
18300 


C 


c: 

o 

c 

c 

2 

c 

c 

c 

c 

c 


INTEGER  N ,  I  ER  ,  I FUL ST ,  I FOL  ,  I  Z  ,  I  V  E  C  ,  N  D  I  V  2  ,  I POI NT 
DOUBLE  PRECISION  MR1 ,mTIlFl,MTILSY,D,Z,wK,SIGN 

DIMENSION  MRU  1  )  ,  MT  I  [_F  L  (  I  FUL  ,  1  )  ,  M  T  I  L  S  V  (  1  )  ,  D(  1  )  ,Z(  IZ,  1)  ,WK(1  ) 


NDI v  2  —  N / 2 


S I  GN  =  1  . 0  DO 


OMPUTE  ORDER  OF  REDUCED  MATRICES 


SET  SIGN  TO  COMPUTE  M~  {  i-  ) 


FIRST  TIME  THRU  LOOP  SOLVE  EIGEN- 
PROBLEM  FOR  M~(+);  SECOND  TIME  FOR 
M-(-) 


10 


2  0 


DO  40  1=1.2 

I POI NT= ( 1-1 ) *NDI V2+ 1 
I F ( I FULST . LT .  t )G0  TO  10 
GO  TO  20 

THEN  COMPUTE  USING  SYMMETRIC  STORAGE  MODE 
CALL  MTLSYMC  NDI V2 . MR  1  , MT I L  S  Y , S I GN ) 

CALL  EIGRStMTILSY , ND I  V  2 ,  IvEC  ,D(  I POI NT )  ,Z(  1  .  I  POINT)  ,  IZ.WK.IER 

) 

GO  TO  30 

ELSE  COMPUTE  USING  FULL  STORAGE  MODE 
I VEC= I VEC+ 1 0 

SIN^.E  IMSL  ROUTINE  E  I  GR  S  REQUIRES 
THAT  MTILFL  BE  EXACTLY  N/2  BY  N/2 

the  fourth  argument  in  the  mtlfl 
call  HAS  BEEN  REPLACED  BY  NDIV2. 

RESTORE  THIS  ARGUMENT  TO  I FUL  IF 
REQUIRED  BY  AN  EIGENVALUE  ROUTINE. 

CALL  MTLFUL(NDIV2 , MR  1  .MTILFl, NDI v2 , SIGN) 

CALL  El GRS (MTI LFL , NDI V2 ,  I VEC  ,  D(  I POI NT )  , Z(  1  ,  I  POINT)  ,  I Z , WK .  I ER 


30 


) 

CONTINUE 


18400 

c 

18500 

c 

18600 

S  I  GN  = 

18700 

c 

1  8800 

40  CONTINUE 

18900 

c 

19000 

RETURN 

19  100 

END 

SET  SIGN  TO  COMPUTE  M~ (- ) 


14 


19  200 
19300 
19400 
19  5  0  0 

19  600 
19700 
198  00 

1  9900 

2  0000 
2  0  10  0 
20200 
203  00 
2  04  00 
20500 
2  060  0 

20  7  00 
208  00 
20900 
2  1  000 
2  1  1  0  0 
2  1  200 
2  1300 
2  1  4  00 
2  15  0  0 
2  16  0  0 
2  17  0  0 


SUBROUTINE  MTLSYM( NDI V2 , MR  1  .  MT  I  LDA  ,  S  I  GN  ) 

Q  ^if***************.************  *******  *  **********  ’************** 

C  ***  THIS  SUBROUTINE  returns  M~(+)  or  m-(-)  stored  IN  *** 

C  ***SYMMETRIC  STORAGE  MODE  IN  THE  ONE-DIMENSIONAL,  DOUBLE  *** 
C  *  * *PREC I S I  ON  ARRAY  MTILDA.  M~(+)  AND  M~  (  -  )  ARE  COMPUTED  *  *  + 
C  *  *  *  FROM  THE  FIRST  ROW  OF  MATRIX  M,  NAMELY  MR  1  ,  WHEN  SIGN  *  *  * 
C  ***IS  +1  AND  -1,  RESPECTIVELY.  *  *  * 

£  +  *  +  +  .i<  +  **  +  :i<  +  *****.i‘* +  +  ****♦*♦  +  *♦ +  *  +  *  +  ********************* 

INTEGER  I  ,  J , K , NDI V2 , N , L 1  , L2 
DOUBLE  PRECISION  MR  1  , MT I LDA . S I GN 
DIMENSION  MR  1  (  1  )  . MTI LDA (  1  ) 

C 

N  =ND I V2*  2 
C 

I  =0 

DO  20  K= 1 , NDIV2 
DO  10  J  = 1  , K 
1  =  1  +  1 
L  1  =  K  -  J  +  1 
L2  =  N  +  2-K~  J 

M TIL DAI  I  ) =  MR  1  ( L 1  )  +  MR  1  ( L2  J  *  SI GN 

10  continue 

20  CONTINUE 

c 

RETURN 

END 


2  1  800 
2  1900 
22000 
2  2  1  00 
2  2  200 
2  2  3  00 
2  2  400 
2  2500 
2  2  6  0  0 
2  2  7  0  0 
2  2  800 
22900 
23000 
2  3  I  00 
2  3  2  0  0 
23300 
23400 
23500 
23600 
23  7  00 
23600 
23900 
24000 
2  4  10  0 
2  4  200 
243  00 


SUBROUTINE  MTLFOL ( ND I  V 2 , MR  1  .MTILDA  ,  I  FUL  , SIGN) 

Q  ************************************************************* 

C  ***  THIS  SUBROUTINE  RETURNS  M~(+)  OR  M-(-)  STORED  IN  *** 

C  *  *  *  F  ULL  STORAGE  MODE  IN  THE  TWO-DIMENSIONAL,  DOUBLE  PRE-  *  *  * 
C  ***CISION  ARRAY  MTILDA.  M~(+)  OR  M~ ( - )  ARE  COMPUTED  FROM*** 
C  *  *  *  THE  FIRST  ROW  OF  MATRIX  M,  NAMELY  MR  1  ,  WHEN  SIGN  IS  *** 

C  ***+1  AND  -1,  RESPECTIVELY.  *** 

Q  *******  +  **  +  *+  **  +  *******:********  +  *  +  *******.1<*******  +  **  +  *  +  +  ***** 

c 

INTEGER  N  ,  NDI V2 ,  I  ,  I FUL , J  .  K  ,  Ll  , L2 
DOUBLE  PRECISION  MR  1  .MTILDA , SI  iN 
DIMENSION  MR1( 1 ) , MTILDA ( I FUL, 1) 

C 

N=NDIV2*2 

c 

DO  2  0  K  = 1  ,  NDI V2 
DO  10  J= 1, K 
L  1  =  K  -  J  +  1 
L2=N+2~K- J 

MTILDA.  (  J  ,  K  )  =MR  1  (Ll  J  +MR  1  (  L  2  )  *  S  I  GN 
MTILDA ( K , J ) =MT I LDA ( J . K) 

10  CONTINUE 
20  CONTINUE 
C 

RETURN 

END 
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